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We use a coarse grained solvent model to study the self assembly of two nano-scale hydrophobic 
particles in water. We show how solvent degrees of freedom are involved in the process. By using 
tools of transition path sampling, we elucidate the reaction coordinates describing the assembly. 
In accord with earlier expectations, we find that fluctuations of the liquid-vapor-like interface sur- 
rounding the solutes are significant, in this case leading to the formation of a vapor tunnel between 
the two solute particles. This tunnel accelerates assembly. While considering this specific model 
system, the approach we use illustrates a methodology that is broadly applicable. 



I. INTRODUCTION 

It is widely accepted that solvent water plays an impor- 
tant role in processes that involve interactions between 
hydrophobic surfaces. For a recent review, see Ref. 
Here, we describe the role solvent plays in the self as- 
sembly of idealized nano-scale hydrophobic monomers, 
specifically, two hard spheres with radii R — lnm. With 
such solutes, the solute surface exposed to solvent is suf- 
ficiently large that the solvation free energy is dominated 
by the solvent-solute surface free energy [j], Q • Accord- 
ingly, the free energy of assembly is roughly 

2irj w (R + r w )[2(R + r w )-d], 

where r w ps 0.14nm is the radius of a water molecule, d is 
the distance between the centers of the two cavities in the 
dimerized state, and j w is the water-vapor surface ten- 
sion. For R — lnm, this estimate gives AE ps — Shk-QT, 
where k# is Boltzmann's constant and T is tempera- 
ture. With such a large binding free energy, there is 
little doubt that two nano-scale hard spheres will form 
long-lived dimers in water. The question we consider 
here is specifically how the dimer forms, what dynamical 
pathways lead to assembly. 

The model we use to address this issue is described 
in the next section. Methods of analyzing its dynamics 
are discussed in Section IlIIl Results and conclusions are 
presented in Section HVl 



II. THE MODEL 

We have adapted the model of ten Woldc and 
Chandler Q. In this model, the solvent is a lattice gas 
through which solute particles are allowed to move con- 
tinuously. Recent workQ has established that the ten 
Wolde-Chandler model does follow from coarse graining 
trajectories of an atomistic model of liquid water. This 
solvent model has two parameters: a lattice grid spac- 
ing, I = 0.21nm, and nearest neighbor coupling constant, 
e = 1.51fceT. With these parameters, the lattice gas has 
the same surface tension and compressibility as water at 
standard conditions. The chemical potential of the lat- 
tice gas we consider is fJ, = /i cocx -|-2.25 x 10 feeT, where 



Mcocx = — 3e is the chemical potential of the lattice gas 
at liquid-vapor phase coexistence. With this choice, the 
model is as close to phase coexistence as is liquid water 
at standard conditions. 

The two solute particles we consider are "ideally" hy- 
drophobic in the sense that their only interactions with 
the solvent is to exclude volume. A solute particle with 
a hard core radius of R will thus exclude the center of 
solvent molecules from occupying a spherical volume of 
radius R + r w . Modest attractions between solute and 
solvent, such as those between oil and water molecules, 
could be added to the model but with little effect on 
solvent fluctuations @ 

The net potential energy for the model can be viewed 
as a free energy that results from integrating out density 
fluctuations occurring on length scales smaller than the 
lattice spacing. For lattice gas plus solutes we take it to 
be 

H[{n k }; {vi}} ps njTtj+^[-^+A/x eg (uj)]nj+?7(r). 

(1) 

where rii = 0, 1 is the lattice-gas variable for cell i. The 
primed summation is over nearest-neighbor pairs of lat- 
tice sites, and U(r) is the pair-potential of interaction 
between the solutes when separated by distance r. The 
variable Vi is the volume of lattice site i occupied by so- 
lute particles, and A/i ex (vi) is the free energy of solvation 
for that volume. It is due to coarse graining that excluded 
volume appears as a soft (rather than hard) constraint 
in the model[6]. Since solvation energy at small length 
scales is proportional to volume, 



Afj, ex (vi) = cvi 



(2) 



where c = 0.6fceT/Z 3 is a constant equal to the excess 
chemical potential per unit volume of a small hard sphere 
in waterjy]. For U(r) we use the WCA potential [3], 



U(r) = 4 e[(cr/r) 12 - (cr/r) 6 ] + e r 
= r 



< 2 1/6 (j(3) 
> 2 1 / 6 ( j(4) 



where a = 2R. 

The system evolves with a stochastic dynamics mov- 
ing the solvent through its discrete configurations and 
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moving the hydrophobic spheres continuously in space. 
The solvent is moved from a time i to a time t + 5ti 
with a Metropolis Monte Carlo [1] sweep. A sweep con- 
sists of M moves, where M is the total number of lattice 
sites in the system. For each move we choose a ran- 
dom lattice site i, and attempt a change of state (i.e. 
rii = — ► rii = 1 or rii = 1 — > rii = 0) with an accep- 
tance probability consistent with the Boltzmann weight 
for the energy function in Eq. [TJ The volumes u, are 
fixed by the positions of the solute particles during this 
step. We associate the sweep time with the physical time 
Sti = 5.0 x 10~ 13 s, as this value approximates the corre- 
lation time for bulk liquid density fluctuation of length 
scale I [i.e., Sti = Z 2 /(47r 2 D), where D is the self diffusion 
constant of liquid water] . Sweeps are performed as either 
collections of changes in state of single lattice sites, or as 
changes in nearest-neighbor pairs of lattice sites to con- 
serve net solvent occupancy. The kinetics depend upon 
which procedure is followed, though the mechanism for 
nano-particle assembly appears unaffected. This point is 
discussed further in Section ITVl 

The hydrophobic particles respond to the forces 



F a ({n}, {n k }) = -V c 



+ f a , (5) 



where the subscript a identifies the particular solute (i.e., 
a =1 or 2) and / is the random force that is the remnant 
of the small length-scale density fluctuations. Because 
small length-scale density fluctuations are to a good ex- 
tent, Gaussian[9j, / should be Gaussian with (/) = 
and (|/| 2 ) = 6kBTj/5t s . We use 7 = 6ttt]R, and 77 is 
the viscosity of liquid water. With these forces, the so- 
lute positions at time t progress to those at time t + St s 
according to 

f a (t + St,) = f a (t) + 6 J±F a ({r ( {t)}, {m}). (6) 
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The value of St s is set as the largest time interval for 
which Eq © maintains detailed balance to an acceptable 
extent. This gives 5t s = 2.8 x lCP 14 s. Trajectories are 
generated by the repeated application of the above proce- 
dure. Crucially, this evolution is reversible and preserve 
the equilibrium distribution consistent with the energy 

To initiate our studies of the model, the solute particles 
were placed randomly inside the system and the solvent 
was allowed to equilibrate around fixed solute particles. 
After the initial solvent equilibration the solute particles 
were allowed to move and trajectories were recorded. The 
simulation cell was a cube with a side length of 6.8nm and 
was periodically replicated in each of the three Cartesian 
directions. 

Free energy surfaces discussed in Section HVl were com- 
puted using umbrella sampling[8] with a bias potential in 
the relevant solvent coordinate and a potential of mean 
force in the relevant particle coordinates. 



III. THEORETICAL METHODS 

Many of the tools we used to analyze the dynamics of 
assembly were borrowed from transition path sampling 
methods. For a review, see Ref. [lOj . As the system 
evolves, the state of the system at time t is given by x t , 
which denotes the collection of variables f\, r*2 and {rii} 
at time t. We define a function h(xt) to be equal to 1 if 
Xt corresponds to a dimerized configuration, and oth- 
erwise. Given an initial configuration, xq (particle posi- 
tions and solvent configuration), we define p(xq, t) as the 
probability that a trajectory initialized from that config- 
uration will be in a dimerized state after an observation 
time r. This probability is analogous to the commitor in 
transition path sampling[10]. While r should not be so 
large as to make p(xq, t) independent of xq, it should be 
large enough to allow for dimerization. This time scale 
separation will be discussed further in Section IIVI The 
collapse probability can be written as, 



p{x, t) 



(S(X - Xg)h(x T )) 
(S(x - Xq)) 



(7) 



where the pointed brackets indicate an equilibrium aver- 
age over initial conditions and is equivalent to an average 
over many trajectories, each of length r. 

The contraction of this probability for the high dimen- 
sional x to that for a lower dimensional q(x) is, 



p(q, 



(p(x ,T)6(q - g(x ))} 
(p(x ,t)) 



(8) 



where q(xo) is the initial value of the coordinate q. The 
g's we have in mind are those that can be good reaction 
coordinates. In practice we compute p(g, r) by averaging 
N trajectories over many sets of initial conditions 



1 N 



(9) 



(i) 

where x T refers to x at time r for the ith trajectory. 

There is an ensemble of transition state configurations, 
and for each member of this transition state ensemble 
p(xq,t) = 1/2. Members of the transition state ensem- 
ble may appear diverse, but there are, however, config- 
urational features that are common to members of the 
transition state ensemble which arise as signatures of the 
reaction mechanism. A reaction coordinate q which prop- 
erly characterizes the dynamics of the reaction should 
also resolve those signature features common to mem- 
bers of the transition state ensemble. If the coordinate 
q does in fact resolve the important features of the re- 
action mechanism, members of the transition state en- 
semble will be narrowly distributed around q — q*, with 
p(q* , r) w 0.5. This expectation, plus the expectation 
that p(xq,t) is also narrowly distributed at q* = q(xo) 
will provide us with a criteria to reject unsatisfactory (or 
poor) reaction coordinates. 
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FIG. 1: Snapshots of a reactive trajectory at several points in time. Each snapshot is labeled on the trajectory plotted in the 
one dimensional phase space of ri2 (top left) and the two dimensional phase space of ri2 and iVvap (bottom left). See text for 
definitions. 



The distribution of p(xq,t) is 



IV. RESULTS AND DISCUSSION 



P(p;q* 



(6{p - p(x , T))S(q* - g(x ))) 
(6(g*-q(x ))) 



where q(xo) is the initial value of the proposed coordi- 
nate. In practice, the distribution is computed by aver- 
aging over many trajectories with a set of N initial con- 
ditions taken from an equilibrium ensemble with q(xo) 
constrained to q* [l2j , 



1 N 

P(p;g*)a-£*(p-p(4V)) 



Figure Q] shows a series of snapshots from a single re- 

(10) 

active trajectory. When viewing trajectories as a func- 
tion of particle separation, r±2 = \r\ — r^j, where f\ and 
r*2 are the positions of the two solute particles, one ob- 
serves three distinct behaviors. At large values of r\%, 
the particles exhibit independent diffusive motion. At 
lower values of r\%, solutes are bound together in a sta- 
ble dimerized state. At intermediate values of 7*12, the 
nano-particles undergo aggregation, which is manifest in 
trajectories as a rapid decrease in 7*12 in time, such a 
rapid decrease is seen in the top left panel of figure Q] 
The time over which this rapid decrease occurs is about 
(11) 50ps in the pictured trajectory. We have found this to 
be typical of all the trajectories we have studied. 
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A. Solute Reaction Coordinate 

Since r%2 is capable of distinguishing between the diffu- 
sive and collapsed regime, it is a natural order parameter. 
To determine if r 12 is also a good reaction coordinate we 
must first define the transition state for assembly. We de- 
fine transition states as those configurations x for which 
p(x, t) = 1/2. This definition requires a specification the 
indicator function h(x) and the time r. 

Judging from the trajectories, we have taken, 

h(x) = 1 ri2 < 2.1nm (12) 
= ri2 > 2.1nm. (13) 

and we take the observation time r = [(7-12 — 
2.1nm)/0.6nm]50ps. The solid horizontal line drawn in 
figure Q] represents the value of 7*12 = r* = 28. 5A for 
which p(r*,r) w 0.5. The top panel of Fig. [2] shows 
P(p; r) for r±2 = T* — 2.85nm. We see that P(p; r) is bi- 
modal, indicating that a given configuration drawn from 
an equilibrium ensemble at 7*12 = r* is typically reac- 
tive (resulting trajectories mostly result in aggregation) 
or nonreactive, but not a reactive thrcshhold. We con- 
clude that r 12 alone does not suffice in describing the 
mechanism of hydrophobic assembly. 

B. Solute and Sovent Reaction Coordinate 

Since the solute particles are spherically symmetric, 
and the solvent is isotropic, the only information con- 
tained in the initial conditions which is not contained in 
7"i2 involves the solvent configuration. The bottom panels 
of Fig. [2]show the average initial solvent density between 
the two solute particles at a separation of r*. The least 
reactive initial solvent configurations have a high solvent 
density in the volume between the two particles while the 
most reactive initial solvent configurations have a vapor 
tunnel which connects the solvent cavities of the two so- 
lute particles. This is reminicent of the vapor tunnel 
which forms as a precursor to the drying of a metastable 
liquid confined between extended hydrophobic surfaces 

[ii mm. 

The difference between panels A and B in Fig. [2] in- 
dicate that a solvent coordinate should be capable of 
resolving a vapor tunnel. A reaction coordinate which 
differentiates between panels A and B and captures the 
interfacial fluctuations is the total number of vapor-like 
(rij = 0) lattice sites in the vicinity of the hydrophobic 
solute particles, A vap , which we define as, 

iVvap = 5^(1 - ni)e(r c - min(|r- - n\, \n - f 2 |)) (14) 

i 

where Q(x) is the Heaviside function (i.e., 1 for x > 
and for x < 0) and r c is a cutoff radius, which we 
have taken to be four lattice spacings beyond the surface 
of a hydrophobic sphere, see Fig. [3J in particular, r c = 
1.84nm. 
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FIG. 2: (top)The distribution of collapse probabilities, P(p; r) 
for ri2 = r* = 2.85nm(lines are guides to the eye), (bot- 
tom) The cylindrical average of the initial solvent configu- 
ration around the solute particle at a distance 7-12 = r* . 
Panel A is an average over the least reactive initial solvent 
configurations (j}(xo, T )\r 12 =r* < 0.10) and panel B is an 
average over the most reactive initial solvent configurations 
(p(xo,r)\r 12 =r* > 0.90). In each panel the center of the par- 
ticles have (rii) = and contour lines are drawn at increasing 
increments of (n») =0.2. 



The bottom left panel of Fig. [T]shows a reactive trajec- 
tory plotted in the two dimensional reaction coordinate 
of 7-12 and A vap along with contours of the free energy 
surface. Focusing first on the contours of the free energy 
surface in Fig. [1] we observe that the diffusive regime is 
shaped like a trough which is parabolic in the direction 
of A vap near the mean (fluctuations in A vap about the 
mean are Gaussian). The diffusive trough in the free en- 
ergy surface funnels into the deep basin associated with 
the aggregated dimer state. The difference in free energy 
between the minimum of the diffusive basin and the bot- 
tom of the aggregated basin is about 64/cbT. Near the 
transition state, fluctuations in A vap are non-Gaussian, 
and large fluctuations, especially towards larger values of 
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N vap are stabilized. In addiction, the assembly of the hy- 
drophobic nano-particles is an activationless process, as 
there is no appreciable barrier in the free energy surface. 
As a result, when the particles are separated by distance 
r* = 2.85nm the solvent fluctuations required to create 
a vapor tunnel are of the order feeT", the thermal energy 
of the system. 

The collapse probability distribution P(p; r, N) was 
computed for two points in the transition state region. 
We found that just as for the one-dimensional ri2 (plot- 
ted in in Fig. \2§ , the distribution of collapse probabilities 
are bimodal for both points in the transition state region 
for the two dimensional coordinate. This is a sign that 
a solvent coordinate with greater spatial resolution is re- 
quired to characterize this transition. Although -/V vap is 
capable of distinguishing between panels A and B in Fig. 
|2l it does not reliably resolve a vapor tunnel in individual 
trajectories. This can be seen by comparing snapshots A 
and C in Fig. [TJ which appear very different in their sol- 
vent configurations, but lie in almost identical regions of 
the two-dimensional phase space (r±2, N mp ). 



C. Further Resolution of Sovent Dynamics 

To achieve greater spatial resolution we restrict our 
observation to the volume between the solute particles, 
where the formation of the vapor tunnel occurs (see 
shaded region in Fig. [3]). The solvent density in this 




■ Observed Volume 

FIG. 3: A two dimensional depiction of the volume (in grey) 
over which the solvent coordinate p cy i is computed. 



region, p cyl is 
where 



Pcyi 



Z^n. t Vj(fi,f 2 ) 



(15) 



9(r c - mindfi - f*i| 5 |r* - r 2 |)) 
xe(min(|r i - fi|, \r t - f 2 |) - R) 
x6(r cyl - \(fi - fi) x fi 2 /ri 2 |) 
x6((r l - r 2 )-ri 2 )9((ri - f$)-fi 2 ) 



with r cy i = 0.7nm. 

Figure [4] plots points visited by many reactive trajec- 
tories along with the contour lines of the free energy sur- 
face for (ri 2 , /?cyl)' With the spatial resolution of p cy \ 




Pcyi 



FIG. 4: Points (black) visited by 100 trajectories projected 
onto the ri2 and p C yi plane, with contours of the free energy 
surface overlaid (white lines). All trajectories were generated 
with initial values of ri2 > 2.8nm. Contour lines are in in- 
crements of 1&bT- Points shown in red are members of the 
transition state ensemble. 

we observe a pronounced depletion of the solvent density 
between the solutes prior to the passing through of the 
transition bottleneck. The red points in figure[2]are mem- 
bers of the transition state ensemble, which are clustered 
near the transition bottleneck of the free energy surface. 



D. Mechanism of Dimerization 

This figure, along with Fig. [2] provides us with a de- 
scription of the aggregation. The transition begins when 
the particles have diffused to positions that bring their 
surrounding solute cavities to within the range accessible 
by equilibrium fluctuation of the cavity interface. As the 
interfaces surrounding each solute particle samples many 
configurations the cavities come into contact which cre- 
ates the beginning of a vapor tunnel. The growth of the 
vapor tunnel is manifest on Fig. 0] as the rapid decrease 
of the density, p C yh between the solutes prior to passing 
through the bottle neck in the free energy surface. The 
existence of a vapor tunnel leads to an unbalancing of 
the solvent induced force on the solute particles. This 
unbalanced force arises because the solute particles are 
being pushed by solvent in every direction accept the di- 
rection of the vapor tunnel. The result is a net solvent 
induced force on each solute in the direction of the vapor 
tunnel. The particles are pushed closer together and are 
eventually squeezed into physical contact by the solvent. 
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E. Role of Solvent Mass Conservation 

Solvent dynamics were carried out along a Markov 
chain, the details of which are provided in Section [Til 
One consequence of using single lattice site moves, is that 
the net solvent density is allowed to fluctuate in a fashion 
that does not account for solvent mass conservation. To 
examine the extent to which this conservation is impor- 
tant, we have generated trajectories with Kawasaki [lCjj 
neighbor exchange dynamics. These dynamics are car- 
ried out just like the Metropolis dynamics described in 
Section [UJ but instead of attempting changes of state 
for single lattice sites, the states of two neighboring lat- 
tice sites are exchanged. The acceptance probability, of 
course, is consistent with the Boltzmann weight for the 
energy in Eq. [JJ The free energy surface, an equilibrium 
property, is invariant to the choice of dynamics. Those 
aspects of trajectories determined by statistics are un- 
changed. Specifically, under both the mass-conserving 
and non mass-conserving dynamics, dimerization is pre- 
empted by the formation of a vapor tunnel which ac- 
celerates assembly. In both cases the onset of acceler- 
ated assembly occurs near ri2 = r* — 2.85nm. But the 
rate at which solutes are pushed together in this region 
are very different for the two forms of solvent dynamics. 
Specifically, after passing through the transition state by 
forming a stable vapor tunnel, trajectories with the mass- 
conserving solvent dynamics take about 2000ps to dimer- 
ize, more than one order of magnitude slower than the 
dynamics without mass conservation. Thus, although 
the mechanism of hydrophobic assembly is invariant to 
solvent mass conservation, the rate of passing through 
the transition state ensemble does depend upon conser- 
vation, so that prefactors to an Arrhenius rate constant 
expression depend upon it significantly. 

F. Role of Solute-Solvent Attractions 

The model considered in detail here contains ideal- 
ized hydrophobic solutes who's only interactions with 



the solvent is to excluded it from occupying a specific 
volume. In real physical systems, there are additional 
solute-solvent interactions, such as Van der Waals attrac- 
tions which can also play a role in the dynamics. Recent 
work[5| has described how the average solvent density 
surrounding hydrophobic solutes responds linearly to at- 
tractive solute-solvent forces. This response increases the 
mean density, and shifts the position of the surrounding 
interface closer to the solute surface. To the extent that 
solvent fluctuations are Gaussian, the nature of the inter- 
facial fluctuations are unchanged. Thus, for fixed solute 
positions, the addition of such attractions necessarily in- 
creases the distance between the interfaces surrounding 
different solute particles. As a result, larger solvent fluc- 
tuations are required to form a vapor tunnel, leading to 
an increase in the activation energy to reach the tran- 
sition state. A free energy barrier to the formation of a 
vapor tunnel that is more than a few k^T , can be reduced 
(and even avoided all together) if the solutes diffuse closer 
together. Thus, attractive forces that do not destroy the 
interface by pinning it to the solute surface (in which case 
the solute could hardly be considered hydrophobic) do 
not qualitatively change the mechanism discussed here. 
The effect of such attractions would, however, cause a 
decrease in the value of the distance r%2 = r* , the dis- 
tance at which the onset of dimerization is observed, in 
effect, small attractive forces make the solvated particles 
look smaller. 
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